function [err, EPSCresult]=EPSC_COV()
% Analyze the Covariance of EPSPs
% Derived from EPSP_TC
% Expects data to be a series of EPSCs in a train, 10 msec apart.
% Analyzes according to Scheuss et al,. 2002.
% Paul B. Manis, Ph.D.
% pmanis@med.unc.edu
% Initial versions: 4/2002, 5/2002.
% 6/2/2002.
%

%try
   sf = getmainselection;
   if(sf > 0) 
      pflag = getplotflag;
      QueMessage('EPSC_COV analysis', 1); % clear the que
      for i = 1:length(sf)
         EPSC_COV2(sf(i), 1, pflag);
      end;
    if(length(sf) >= 3)
        EPSC_mean_var;
    end;
   end;
   return;
   %catch
   watchoff;
   QueMessage('Error in EPSC_COV Analysis routine', 1);
%
%-----------------------------------------------------------

function [err, EPSCresult]=EPSC_COV2(sf, plot_flag, pflag)

global VOLTAGE CURRENT DFILE ALLCH
global CONTROL
% first generate empty output arrays

err = 0;
EPSCresult=[];

dat = [];
time = [];
ivresult = []; % initialize it
do_spike = 0;	% 0 turns off detailed spike analysis; 1 turns it on

extrap = 1; % extrapolate from previous response to next response... to remove the baseline.

% this next routine reads all of the data and information we need,
% including calculating some standard stuff that is put into DPAR.
% you will want DPAR.time especially - it is the time base that goes
% with the data.

[DFILE, DPAR, err] = ana_setup(DFILE, sf);
if(err ~= 0)
    return;
end;
isrc = 1;

if(~isempty(ALLCH))
    VOLTAGE = ALLCH{isrc};
    CURRENT = ALLCH{isrc+1};
end;



% Get analysis windows 
QueMessage('EPSC_COV - analysis windows');
stim_list=CONTROL(sf).stim_time; % get the array
psp_time=CONTROL(sf).psp_time; % get psp definition array
if(ischar(psp_time))
	psp_time = str2num(psp_time);
end;
if(ischar(stim_list))
	stim_list= str2num(stim_list);
end;
%l = length(psp_time);
%if(l > 4)
%	psp_time = reshape(psp_time, l/2, 2)';
%end;
p2=size(psp_time);
npsc = length(stim_list); % get number of pscs to examine

trmp = floor(number_arg(stim_list(1))./DPAR.RATES - 0.5); % for RMP/ihold determination
psp_time = reshape([stim_list+0.5, stim_list+4],length(stim_list),2);

for i = 1:npsc % compute time windows for each psc in the train
   t0(i,:) = floor(number_arg(psp_time(i, 1))./DPAR.RATES);
   t0b(i,:) = floor(number_arg(psp_time(i, 1)-0.2)./DPAR.RATES);
   t1(i,:) = floor(number_arg(psp_time(i, 2))./DPAR.RATES);
   st(i,:) = floor(number_arg(stim_list(i))./DPAR.RATES);
end;
   
%-----------------------------Begin analysis in earnest --------------------------------
QueMessage('EPSC_COV - Measuring Ih and rmp');

% measure holding current and "rmp"
% also measure Rin with s1 pulse after ho (ts2)
% and tau by fitting to s1 pulse.
for i = 1:DPAR.records
	hold_cur(i) = mean(CURRENT(i,1:trmp(i))); % save in array for later usage
	RMP(i) = mean(VOLTAGE(i,1:trmp(i)));
end;
CONTROL(sf).iHold = mean(hold_cur);
CONTROL(sf).Rmp = mean(RMP);
[CURRENT] = FP_artsupp(CURRENT, DFILE, sf);
% Now smooth the Current out a bit
QueMessage('EPSC_COV - Smoothing');
for i = 1:DPAR.records
	fsamp = 1000/DPAR.RATES(i); % get sampling frequency
	fco = 10000;		% cutoff frequency in Hz
	wco = fco/(fsamp/2); % wco of 1 is for half of the sample rate, so set it like this...
	if(wco < 1) % if wco is > 1 then this is not a filter!
	   [b, a] = fir_win(8, wco); % fir type filter... seems to work best, with highest order min distortion of dv/dt...
	   ismo(i,:) = DigitalFilt(b, a, CURRENT(i,:)')'; % filter all the traces...
	else
	   ismo(i,:) = CURRENT(i,:);
	end
end

QueMessage('EPSC_COV - Finding EPSC peak currents');% find EPSC amplitudes for each entry in the train
absipscpeak = 0;

for r = 1:DPAR.records
    %   ismo(r,:) = ismo(r,:) - ibase(r);
    ismo(r,:) = ismo(r,:) - hold_cur(r);
    twofiftyusec = floor(0.250/DPAR.RATES(r)); % number of points in 250 usec...
    lastjmin = 0;
    for i = 1:npsc
        [imin, jmin] = min(ismo(r,t0(i,r):t1(i,r)));
        imin = mean(ismo(r, (t0(i,r)+jmin-3):(t0(i, r)+jmin+3)));
        if(imin < absipscpeak)
            absipscpeak = imin;
        end;
        if(i > 1 && i <= npsc && extrap)
            % use chebyshev to get a fit to the trace from the peak time +0.5 msec to the next stim -0.5 msec
            ts = [0:fitwinr-fitwinl]*DPAR.RATES(r);
            fitdata = ismo(r, fitwinl:fitwinr);
            [a0 a1 tau] = expfit(ts, fitdata);
            [a b] = linreg(ts, fitdata);
            % now calculate for our current min
            tmin = (stim_list(i)-stim_list(i-1)) + jmin*DPAR.RATES(r);
            iextrap = a0 + a1*exp(-tmin/tau);
            iextrap2 = a*tmin+b;
            tf = a0 + a1*exp(-ts/tau);
            tf2 = a*ts+b;
            htp = newfigure('testplot', 'testplot');
            plot(fitdata);
            hold on;
            plot(tf, 'r-');
            plot(tf2, 'g-');
            err1 = (fitdata-tf).^2;
            err2 = (fitdata-tf2).^2;
            if(err1 < err2)
                I_psc(i, r) = imin - iextrap;
            else
                I_psc(i, r) = imin - iextrap2;
            end;
                
        else
            I_psc(i, r) = imin; %  - mean(ismo(r, t0b(i, r):t0(i,r)));
        end;
        fitwinl = t0(i,r)+jmin+twofiftyusec;
        if(i < npsc)
            fitwinr = st(i+1,r)-twofiftyusec;
        else
            fitwinr = st(i)+(st(i)-st(i-1))-twofiftyusec;
        end;
    end;
end;

% reduce to a list without failures - must be true for all pulses in consequective records.
r_ok = [];
failthresh = number_arg(CONTROL(sf).thresh);
for r = 1:DPAR.records-1
   if(all(abs(I_psc(:, r)) > failthresh) && all(abs(I_psc(:, r+1)) > failthresh))
      r_ok = [r_ok r];
   end;
end;



% compute mean current for each stimulus pulse
for i = 1:npsc
   I_m(i) = mean(I_psc(i,r_ok));
   I_var(i) = var(I_psc(i,r_ok));
end
ipscmin = min(min(I_psc));

% and compute the variance and the covariance of the amplitudes for each pair...
vsum = zeros(npsc, 1);
covsum = zeros(npsc, 1);
covar = zeros(npsc, DPAR.records-1);
var1 = zeros(npsc, 1);
var2 = zeros(npsc, 1);
j = 0;
if(isempty(r_ok))
    fprintf(1, 'EPSC_COV: r_ok is zero at about line 155? \n');
    return;
end;
for i = 1:npsc
   covsum(i) = 0;
   for r = r_ok
   if(i < npsc)
      j = i+1;
   else
      j = i-1;
   end;
      vsum(i) = vsum(i) + 0.5*((I_psc(i,r)-I_psc(i,r+1))^2);
      covar(i,r) = 0.5 * (I_psc(i,r)-I_psc(i,r+1))*(I_psc(j,r)-I_psc(j,r+1));
      covsum(i) = covsum(i) + covar(i,r);
   end;
   var1(i) = mean(covar(i,:).^2)-mean(covar(i,:))^2;
   var2(i) = (3.5/DPAR.records)*(std(I_psc(i,:))^2)*(std(I_psc(j,:))^2);
end;
rnum = length(r_ok)-1;
vsum = vsum ./ rnum;
covsum = covsum ./ rnum; % mean covariance

Cvq = 0.32; % quantal variability; is CV of mEPSCs - use to correct per Scheuss, eq. 7.
for i = 1:npsc
   switch i
   case 1
      q(i) = (1+Cvq)^2*((vsum(i)/I_m(i)) - (covsum(i)/I_m(i+1)));
	case npsc
      q(i) = (1+Cvq)^2*((vsum(i)/I_m(i)) - (covsum(i-1)/I_m(i)));
   otherwise
      a = (1+Cvq)^2*((vsum(i)/I_m(i)) - (covsum(i)/I_m(i+1)));
      b = (1+Cvq)^2*((vsum(i)/I_m(i)) - (covsum(i-1)/I_m(i)));
      q(i) = (a+b)/2;
   end;
end;
Ncov = - floor((I_m(1)*I_m(2))/covsum(1)); % N is only calculated for the FIRST pair of responses.
Qc = floor(I_m./q); % This is quantal Content for each response

EPSC_Cov.Qc = Qc;
EPSC_Cov.Ncov = Ncov;
EPSC_Cov.Im = I_m;
EPSC_Cov.Ivar = I_var;
EPSC_Cov.npsc = npsc;

CONTROL(sf).EPSC_Cov = EPSC_Cov; % copy the structure over

QueMessage('EPSC_COV - analysis complete');


%-----------------------------Prepare for plotting--------------------------------
% for plotting, do baselines/stddev.
%----------------------------- plot if figure is set--------------------------------
if(plot_flag >= 0)
   hepsc = findobj('Tag', 'EPSC_COV'); % check for pre-existing window
   if(isempty(hepsc)) % if none, make one
      hepsc = figure('Tag', 'EPSC_COV', 'Name', 'EPSC Covariance Analysis', 'NumberTitle', 'off');
   end
   figure(hepsc); % otherwise, select it
   clf; % always clear the window...
   fsize = 7;
   msize = 3;
   tmax = max(time);    
   
   [m, b, r, p] = linreg(r_ok', -I_psc(1,r_ok)'); % check drift in response.
   yline = m*[1:DPAR.records]+b;
   drift = 100*((m*DPAR.records+b)-(m+b))/abs(m+b);

% plot Max voltage for EPSP
   subplot('Position',[0.07,0.35,0.45,0.3]);
   sym = {'rs', 'bo', 'g^', 'k+', 'cd'};
   for i = 1:npsc
      plot(r_ok, -I_psc(i,r_ok)/1000, sym{i}, 'Markersize', 2, 'Markerfacecolor', 'k') % data in black
      hold on;
      if(i == 1)
         plot([1:DPAR.records], yline/1000, 'linestyle', '-', 'color', 'r');
      end;
   end;
   
	set(gca, 'FontSize', fsize);
   ylabel('EPSP Amplitude (nA)', 'FontSize', fsize);
   set(gca, 'XTickLabelMode', 'Manual');
   %grid;
   u=get(gca, 'YLim');
   % plot the EPSPs...
   subplot('Position', [0.07, 0.7, 0.45, 0.25]);
   plot(DPAR.time(r_ok,:)', ismo(r_ok,:)'/1000);
   set(gca, 'XLim', [0, 60]);
   v = absipscpeak*1.1/1000;
   set(gca, 'YLim', [v -v/4]);
   xlabel('ms', 'FontSize', fsize);
   ylabel('nA', 'FontSize', fsize);
	set(gca, 'FontSize', fsize);
   
   % plot 1st vs second epsc
   subplot('Position', [0.57, 0.7, 0.18, 0.25]);
   plot(-I_psc(1,r_ok), -I_psc(2,r_ok), 'ks', 'Markersize', 2, 'Markerfacecolor', 'k', 'linestyle', 'none');
   xl = get(gca, 'XLim');
   yl = get(gca, 'YLim');
   m1 = mean(-I_psc(1,r_ok));
   m2 = mean(-I_psc(2,r_ok));
   hold on
   plot(xl, [m2, m2], 'linestyle', '--');
   plot([m1, m1], yl, 'linestyle', '--');
   xlabel('EPSC1', 'FontSize', fsize);
   ylabel('EPSC2', 'FontSize', fsize);
   set(gca, 'FontSize', fsize);
   
   subplot('Position', [0.8, 0.7, 0.18, 0.25]);
   hp = errorbar([1:npsc], -I_m, sqrt(I_var), 'ko');
   set(hp, 'markersize', 2, 'markerfacecolor', 'k');
   xlabel('EPSC #', 'FontSize', fsize);
   ylabel('EPSC nA', 'FontSize', fsize);
   set(gca, 'FontSize', fsize);
   if(length(I_m) > 1)
       P = I_m(2)/I_m(1);
   else
   P = 0.0;
   end;
   if(length(I_m) >=5)
       P5 = I_m(5)/I_m(1);
   else
       P5 = 0.0;
   end;
   
   % compute parabolic fit to mean-variance plot as follows:
   % 1. must go through (0,0), therefore,
   %    y = Ax^2+Bx; C = 0.
   % 2. transform to linear:
   %    y/x = Ax + B. (CV vs mean)
   % 3. do linear fit.
   
   [A, B, r, p] =linreg(-I_m/10^3, -I_var./(10^6*I_m/10^3));
	   
   % plot of mean vs variance for each EPSP
   subplot('Position', [0.57, 0.4, 0.18, 0.18]);
   plot(-I_m/10^3, I_var/10^6, 'ks', 'Markersize', 2, 'Markerfacecolor', 'k', 'linestyle', '-');
   xl = get(gca, 'XLim');
   yl = get(gca, 'YLim');
   hold on
   xp = [0:0.05:xl(2)];
   plot(xp, A*xp.^2+B*xp, 'k-');
  
   ylabel('Var (nA^2)', 'FontSize', fsize);
   xlabel('Mean (nA)', 'FontSize', fsize);
   xl = get(gca, 'XLim');
   yl = get(gca, 'YLim');
   set(gca, 'XLim', [0 xl(2)]);
   set(gca, 'YLim', [0 yl(2)]);
   set(gca, 'FontSize', fsize);
   
   % plot of mean vs variance/mean for each EPSP
   subplot('Position', [0.80, 0.4, 0.18, 0.18]);
   plot(-I_m/10^3, -I_var./(10^6*I_m/10^3), 'ks', 'Markersize', 2, 'Markerfacecolor', 'k', 'linestyle', '-');
   ylabel('Var/Mean (nA)', 'FontSize', fsize);
   xlabel('Mean (nA)', 'FontSize', fsize);
   xl = get(gca, 'XLim');
   yl = get(gca, 'YLim');
   set(gca, 'XLim', [0 xl(2)]);
   set(gca, 'YLim', [0 yl(2)]);
   set(gca, 'FontSize', fsize);
 
   % plot summary of covariances for pairs...
   subplot('Position', [0.07, 0.07, 0.2, 0.2]);
   plot([1:npsc-1], covsum(1:npsc-1)/10^6, 'ks', 'Markersize', 2, 'Markerfacecolor', 'k', 'linestyle', '-');
   ylabel('Cov (nA^2)', 'FontSize', fsize);
   xlabel('EPSC Pair', 'FontSize', fsize);
  	set(gca, 'FontSize', fsize);
 
   subplot('Position', [0.3, 0.07, 0.2, 0.2]);
   plot([1:npsc], -q, 'ks',  'Markersize', 2, 'Markerfacecolor', 'k', 'linestyle', '-');
   yl = get(gca, 'YLim');
   set(gca, 'YLim', [0, yl(2)]);
   ylabel('quantal size', 'FontSize', fsize);
   xlabel('EPSC #', 'FontSize', fsize);
  	set(gca, 'FontSize', fsize);
 
   subplot('Position', [0.53, 0.07, 0.2, 0.2]);
   plot([1:npsc], Qc, 'ks',  'Markersize', 2, 'Markerfacecolor', 'k', 'linestyle', '-')
   ylabel('quantal content', 'FontSize', fsize);
   xlabel('EPSC #', 'FontSize', fsize);
   set(gca, 'FontSize', fsize);

   % plot textual information abstracted from analysis...
   subplot('Position',[0.74,0.07,0.25,0.25])
  	axis([0,1,0,1])
   axis('off')
   text(0,0.9,sprintf('%-12s R[%d:%d] Protocol:%-8s',DFILE.filename, DFILE.frec, DFILE.lrec, CONTROL(sf).protocol), 'Fontsize', 8);
   text(0, 0.8, sprintf('Comment: %s', DFILE.comment), 'Fontsize', 8);
   text(0,0.7,sprintf('Sol:%-12s  Gain:%4.1f LPF:%4.1f kHz', CONTROL(sf).solution, DFILE.igain, DFILE.low_pass(1)), 'FontSize', 8);
	text(0, 0.6, sprintf('Drift: %6.2f%%', drift), 'Fontsize', 8);
   text(0, 0.5, sprintf('quantal size: %6.1f %6.1f %6.1f %6.1f %6.1f', q(1), q(2), q(3), q(4), q(5)), 'Fontsize', 8);
   text(0, 0.4, sprintf('Qc: %d, %d, %d, %d, %d, Sum: %d, NCov: %d', Qc(1), Qc(2), Qc(3), Qc(4), Qc(5), sum(Qc), Ncov), 'Fontsize', 8);
   text(0, 0.3, sprintf('Facil (p2/p1): %6.3f  (p5/p1): %6.3f', P, P5), 'Fontsize', 8);
   text(0, 0.2, sprintf('EPSC1 mean: %6.3f  var: %6.4f nA^2', I_m(1)/1000, I_var(1)/10^6), 'Fontsize', 8);
   orient landscape
   drawnow
   % control printing and closing of window for automatic runs
   % f = 1 creates plot and leaves it up
   % f = 2 creates plot but closes it when done
   % f = 3 creates plot and prints it and then closes it
   if (pflag > 0)
        eval (sprintf('print -f%d', findobj('tag', 'EPSC_COV')));
   end
   if pflag == 2  
      close(hepsc);
   end
end;

